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, Abstract The elfect of interfacial slip on steady-state and time-periodic flows 

of monatomic liquids is investigated using non-equilibrium molecular dynamics 
simulations. The fluid phase is confined between atomically smooth rigid walls 
and the fluid flows are induced by moving one of the walls. In steady shear flows, 
the slip length increases almost linearly with shear rate. We found that the velocity 
profiles in oscillatory flows are well described by the Stokes flow solution with the 
slip length that depends on the local shear rate. Interestingly, the rate dependence 
of the slip length obtained in steady shear flows is recovered when the slip length 
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1 in oscillatory flows is plotted as a function of the local shear rate magnitude. For 

" both types of flows, the friction coefficient at the liquid-solid interface correlates 
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well with the structure of the first fluid layer near the solid wall. 
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1 Introduction 



The rational design of micro- and nanofluidics devices requires an accurate pre dic 



20051) . It 



tion of time-dependent flows at the submicron scales (jKarniadakis et al. 
is well recognized that fluid flows in confined systems can be significantly affected 
by slip boundary conditions. The velocity discontinuity is usually quantified via 
the slip length, which is defined as an extrapolated distance with respect to the 
liquid-solid interface to the point where the relative tangential velocity component 
vanishes. At sufficiently high oscillation frequencies, the fluid slip velocity might 
not be in phase with t he substrate velocity, and, therefore, the slip len gth in general 



is a complex number (jWillmott and Tallon 



2007 



Ng and Wang 



20111) . Oscillatory 



flows with slip boundary conditions of Ne wtonian liquids were studied experimen- 



tally 



2003; 



using quartz cr ystal microbalance (jFerrante et al 



Du et al 



1994; 



Ellis and Havward 



20041) . Molecular dynamics (MD) simulations are particularly well 



suited to investigate the e ffects of materials properties of liquid-solid interfaces on 



flow boundary conditions (jBocauet and Barrat 



2007 



Li et al 



2010) 



A direct comparison between M P simulations and continuum analysis of steady- 



state flows over chemically textured (Pri eziev et al 



2011 ) or periodically corrugated (jPrieziev and Troia n 2006 



2005 


Oian et al. 


2005 


Prieziev 



Niavarani and Prieziev 



2008a) surfaces has indicated that there is an excellent agreement of the velocity 
profiles and the effective slip lengths when the typical length scale of substrate 
inhomogeneities is about an order of magnitude larger than the molecular size. In 
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the case of time-dependent flows, the main difficulty in extracting flow properties 
using MD simulations is that averaging over thermal fluctuations has to be re- 
peated over many cycles, which often requires significant computational recourses. 
The MD simulations of oscillatory flows of monatomic and polymeric fluids have 
shown that the velocity profiles with no-slip boundary co nditions can be well de- 



scribed by the continuum mechanics (jKhare et al 



20011) . The flow profiles with 



a finite slip velocity in oscil latory flows were reported at low fluid densit ies and 



weak wall-fluid interactions (Han sen and Ottesen 



2006 



Hansen et al 



2007 ). More 



recently, it was shown that the slip length depends on the magnitude and gra- 
dient of shear rate near the oscillating wall, a nd the fluid slip velocity lags th e 



2012) 



wall velocity, thus leading to a hysteresis loop (jThalakkottor and Mohsenil k 
However, the dependence of the slip length on the local shear rate or oscillation 
frequency and amplitude has not yet been systematically investigated. 

Molecular dynamics simulations by iThompson and Troianl (| 1997T) have shown 



that in steady shear flow of Netwonian liquids over atomically smooth crystalline 
surfaces, the slip length is constant at relatively low shear rates and it increases 



slip length was repeatedly observe 


i in MD studies ( 


Prieziev and Troian 


2004 




Yane and Fang 


2005; 


Prieziev 


2007a 


Asproulis and Drikakisi2010i 


Niavarani and Prieziev 


2010; 


Wang and Zha 


3 2011a: 


Pahlavan and Freunc 


2011; 


Wane and Zhao 


2011b 




Kannam et al. 


2012; 


Prieziev 


20121) 


It was also found that the slip length varies 



almost linearly with shear rate when liquid and solid phases form incommensu- 
rable structures at th e interface and the wall-fluid interaction energy is sufficiently 



high (jPrieziev 



2007allbT) . More recently, it was demonstrated that the characteristic 



slip velocity associated with the onset of the nonlinear slip regime correlates well 
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with the diffusion time of fluid monomers over the distance between the nearest 



20101) . One of the 



minima of the periodic surface potential at equilibrium (|Prieziev 
motivations of the present study is to examine whether shear rate dependent slip 
boundary conditions observed in steady-state flows are valid for time-dependent 
flows. 

For steady-state flows, a number of previous MD studies have established a 
correlation between th e degree of slip and fluid structure induced by the perl 



Prieziev 



2007a; 



odic surface potential ([T hompson and Robbins 



Niavarani and Prieziev 



2008b 



1990; 



Barrat and Bocquet 



Prieziev 



2009 



2010; 



1999; 



Zhang et al 



20121) . In particular, it was shown that for atomically smooth, weakly attractive 
surfaces, the friction coefficient at the liquid-solid interface is well described by a 
function of the product of the main peak in the static structure fact or and the 



20101) . However, 



contact density, both evaluated in the first fluid layer ([Priezievl 
the situation is less clear for time-dependent flows where the surface-induced fluid 
structure and boundary slip might have a phase difference (especially at high oscil- 
lation frequencies), and the conclusions obtained for steady-state flows might not 
be valid. In the present study, we performed a comparative analysis of the fluid 
structure and the friction coefficient for steady-state and time-periodic flows. 

In this paper, non-equilibrium steady-state and time-periodic molecular dy- 
namics simulations are performed to investigate Couette flows with slip boundary 
conditions. First, the rate dependence of the slip length is computed in steady 
shear flows. Then, the velocity profiles in oscillatory flows are compared with the 
Stokes flow solutions in a wide range of frequencies. We find that the slip length as 
a function of the local shear rate estimated at the stationary and oscillating walls 
is in good agreement with the results obtained for steady flows. We will also show 
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that, for both types of flows, the friction coefficient at the liquid-solid interface 
correlates well with the structure of the fluid layer in contact with the solid wall. 

The rest of the paper proceeds as follows. In the next section, the details of 
molecular dynamics simulations are described. In Sect. O we briefly review the 
Stokes flow solution for oscillatory Couette flows, present the results for steady- 
state shear flows, and then analyze the velocity and density profiles, slip length, 
and fluid structure in oscillatory flows. Conclusions are given in the last section. 

2 Molecular dynamics simulation model 

The model system consists of Nf = 4608 fluid monomers confined between rigid 
atomistic walls as shown in Fig. [I] The pairwise interaction between any two fluid 
monomers is modeled by the Lennard-Jones (LJ) potential 



where e and a are the energy and length scales, and the cutoff radius r c = 2.5 a. In 
our simulations, the same parameters are used to describe the interaction between 
fluid monomers and wall atoms; namely, e w f = e, cr w f = a, and r c = 2.5 a. The 
wall atoms are fixed rigidly at the lattice sites and do not interact with each other. 

The viscous heating generated in the oscillating flow was removed by means 
of the Langevin thermostat, which was coupled only to the equation of motion 
perpendicular to the plane of shear as follows: 




(1) 




(2) 



(3) 




(4) 
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where r = 1.0 r 1 is the friction coefficient and h is a random force with zero m ean 



1990). In 



and variance (fi(0)fj{t)) = 2mkBTLrS(t)5ij (jThompson and Robbing 
our setup, the temperature of the Langevin thermostat is set Tl = 1.1 s/ks, where 
kg is the Boltzmann constant. The equations of motion wer e integrated numeri- 
cally using the fifth-order Gear predictor-corrector algorithm (jAllen and Tildeslev 



19871) with a time step At = 0.005 r, where r = \Jma 2 js is the characteristic 
LJ time. The length, energy, and tim e scales for liquid argon a re a = 0.34 nm, 



e/k B = 120 K, and r = 2.16 x 10" 12 s (jAllen and Tildeslev 



1987) 



The fluid phase of density p = 0.81 cr -3 is confined between two crystalline 
walls with density p w = 2.73 cr -3 , as illustrated in Fig.[T] Each wall consists of two 
layers of atoms arranged rigidly on sites of the face-centered cubic (fee) lattice. The 
lateral dimensions in the xy plane are measured L x = 25.03 a and L y = 9.63 a, 
and the channel width is fixed h = 23.58 a. Periodic boundary conditions were 
applied in the xy plane parallel to the solid walls. This simul ation setup is v ery 



2007alM 



similar to the one used previously for steady Poiseuille flows (|Priezievl : 
except that in the present study the system size in the y direction is slightly larger. 
In steady shear flows, the fluid viscosity p, = (2.15 ± 0.15) era~ 3 was found to be 
shear rate and temperature in dependent for 7T < 0.072 and 1.1 ^ Tks/e ^ 1-35 



(jNiavarani and Prieziev 



2010) 



To simulate the oscillatory Couette flow, the upper wall velocity was varied in 
the x direction with the angular frequency ui and amplitude U, while the lower wall 
always remained stationary. In the present study, the oscillation frequency was set 
ujt = 10 _1 , 10 -2 , 10 -3 , and 10 -4 (see Table[l|. Before the averaging procedure, 
the steady-periodic flow was equilibrated during the time interval of about 5 x 



10 t. The measurements of the velocity, density, and temperature profiles were 
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made at discrete times cjt n = mr/A + 2nm, where n = 0, 1, 7 and m is the 
integer. These profiles were averaged within horizontal bins of thickness Az = 
0.01 a during the time interval T/100, where T = 2-k/uj is the period of oscillation. 
A typical simulation time at low shear rates is about 1.2 x 10 6 r. 



3 Results 



3.1 Hydrodynamic predictions 



The problem of fully-developed oscillatory viscous flow confined betw een two paral- 



l el wal ls wit h slip boundary con d itions was considered analytically by 



Khaled and Vafai 



(|2004f ) and 



Matthews and Hill 



(2009). Below, we briefly review the problem and 
its solution for the flow geometry depicted in Fig.[T] The x-component of the mo- 
mentum equation (parallel to the walls) is given by 



du x 

'"aT 



d 2 u x 



(5) 



where fi and p are the fluid viscosity and density. The boundary conditions at the 
top and bottom walls are specified as follows: 



ztop ■ u x = Usm(wt) - 

OZ 



Z = Zbot : U x = L; 



du x 



(6) 
(7) 



where U is the amplitude and u is the frequency of oscillation. The slip lengths at 



the u pper and lower walls L\ ^ Li are assumed to be constant (Matthews and Hill 



2009). We note that the special case L\ = L2 was considered by 



Khaled and Vafai 



(2004). 
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The solution of the problem Eq. (0 subject to the boundary conditions Eqs. ([6]- 
[7]) is given by 

u x {£) = A2 U +B2 (exp(+Kz) x { [L 2 K{A + B) + A] sin(urt + K z) 

+ [L 2 K(A -B)-B] cos(wi + Kz)} 
+exp(-Kz) x { [L 2 K(A + B) - A) sin(wf - Kz) 

+ [L 2 K{A-B)+B]coa(ujt-Kz)}y (8) 

where 

A = A + — A~ and B = B + + B~ , (9) 
which in turn are defined as follows: 

A ± = exp(±Kh) { [l ± (Li + L 2 )K] cos(Kh) 

-[(Li + L 2 )K ±2L X L 2 K 2 ] sm(Kh)}, (10) 

B ± = exp(±Kh) { [l ± (Li + L 2 )K] sm(Kh) 

+ [(Li + L 2 )K ±2LxL 2 K 2 ] cos(Kh)}, (11) 

and K = ^ujp/2fj,. In Sect. 13.31 the velocity profiles obtained from MD simulations 
will be fitted to Eq. (JSJ) with the parameters L\ and L 2 . 

3.2 Steady shear flows 

The simulations were first performed at steady-state flow conditions when the 
upper wall was translated with a constant velocity, while the lower wall always 
remained stationary. The upper wall velocity was varied in the range 0.025 tr/r ^ 
U ^ 6.5 a jr. The lower limit was chosen to reduce the averaging time due to 
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thermal fluctuations, while the upper limit was set to avoid the nonlinear slip 
regime at very high shear rates when the slip velocity becomes much larger than 



the fluid thermal velocity (jNiavarani and Prieziev 



20101 ) . In the present study, the 



maximum slip velocity and shear rate in steady shear flows are about 1. 59 cr/r 
and 0.14t -1 , respectively, which provide an upper estimate of the Reynolds num- 
ber Re ~ 29.5. It was previously shown for slip flows over periodically corru- 
gated surfaces that the inertia term in the Navier-Stokes equation produces a 
noticeable difference in the slip length at higher Reynolds numbers of about 130 



(jNiavarani and Prieziev 



2008a) 



The representative velocity and density profiles for the upper wall speeds 
U = 0.1 a /t and U = 4.0 <t/t are plotted in Fig. [2] As expected, the fluid den- 
sity profiles exhibit pronounced oscillations near solid walls that gradually decay 
to the uniform bulk value. The magnitude of the first density peak defines the 
contact density p c . Notice that the amplitude of the density oscillations is slightly 
reduced at the higher upper wall speed U = 4.0 cr/r. The corresponding velocity 
profiles, normalized by the upper wall speed, are linear throughout the channel 
and are characterized by the finite slip velocity at both walls. As clearly observed 
in Fig.[2](b), the relative slip velocity is larger at the higher upper wall speed. Also, 
it was shown previously for a similar MD setup, that the fluid temperatu re near 



2007allbf) . The 



the interfaces increases by about 10% at high shear rates (jPriezievl l' 
correlation between the contact density and fluid temperature in the first layer as 
a function o f the slip veloc ity was recently reported for polymeric fluids in steady 



shear flows (jPrieziev 



2012) 



For steady-state flows, the slip length was estimated from the linear extrapola- 
tion of the the velocity profiles to u x (z) = below the lower wall and to u x (z) = U 
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above the upper wall, and then the two values were averaged. The variation of the 
slip length as a function of shear rate is presented in the inset of Fig. 0(a). In agree- 
ment with the results of previous studies, where the behavior of the slip length 
was investigated in a wide range of shear rates and wall-fluid interaction energies 



(jPrieziev 



2007aUbh , the slip length increases almost linearly with shear rate when 



e w f = s. It is expected, however, that when the wall-fluid interaction energy is 



reduced (jPrieziev 



2007aTI . then the magnitude of the slip length increases and its 



shear rate dependence can b e well fitted by the power-law function proposed by 



Thompson and Troiarj 



1997). In the next section, the boundary conditions and 
fluid structure computed in steady shear flows will be compared with the results 
obtained for time-periodic flows. 



3.3 Oscillatory Couette flows 

We next consider oscillatory flows driven by the upper wall, u™(t) = {7sin(w£), 
with frequencies lot = 10 -1 , 10 -2 , 10 -3 , and 10 -4 . The corresponding period and 
amplitude of oscillations, the Stokes boundary layer thickness, as well as the upper 
estimate of the Reynolds numbers are given in Table[T] For each frequency, the 
amplitude of the velocity oscillation U was chosen such that the fluid slip velocity 
at the upper wall was always less than 1.5 a jr. As can be seen from Tabled 
the thickness of the Stokes boundary layer is smaller than the channel width at 
higher frequencies lot = 10 -1 and 10 -2 . Nevertheless, the upper estimate of the 
Reynolds number, based either on the channel width or the Stokes layer thickness, 
is about 27.6, which is indicative of laminar flow conditions. In the present study, 
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the smallest amplitude of the upper wall velocity was set U = 0.25 a jr in order to 
compute accurately the velocity profiles without excessive computational efforts. 

Examples of the velocity profiles for different frequencies lo and amplitudes U 
are presented in Figs.O |4] [5] [6] and [7] The MD data were averaged over about 20 
periods at the lowest frequency lot = 10 -4 and over 2 x 10 4 periods at the highest 
frequency lot = 10 _1 . In all figures, the red curves represent the least square fits 
of the MD data to Eq. (JSj) with the parameters Li and L2- The shear rate at the 
upper and lower walls was then computed by taking the derivative of the best fit 
function du x /dz at z = ±11.79 a. We found that the MD velocity profiles are well 
described by the continuum solution Eq. ([8]), except in the interfacial regions of 
about 2 a at high shear rates. The small discrepancy observed between the MD 
and continuum results may originate from the inertial effects and/or the fluid 
temperature increase near the walls at high shear rates. 

At the highest frequency lot = 0.1, the Stokes layer thickness is nearly three 
times smaller than the channel width (see Table[l|, and, therefore, the velocity 
profiles near the stationary lower wall are not significantly affected by the moving 
upper wall, and, as a result, the interfacial shear rate at the lower wall remains 
relatively low (see Fig. [3}. When lot = 0.01 in Fig. 21 the Stokes layer thickness 
is approximately equal to the channel width, and the slip velocities and shear 
rates at the upper and lower walls become comparable. Furthermore, at lower 
frequencies, lot = 10 -3 and 10 -4 , the flows appear to be quasi-steady and the 
velocity profiles are nearly linear throughout the channel (see Figs. [6] and [7]). At the 
lowest frequency lot = 10 -4 , the velocity profiles are almost indistinguishable when 
the magnitude of the upper wall velocity is the same (e.g., when cot = 7r/4 and 3 7r/4 
in Fig.[JJ). Notice also that in all cases except lot = 10 -4 , the velocity profiles at 
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times tot = and 7r are symmetrical to each other with respect to the line u x = 
and the fluid slip velocity is not zero, indicating that the fluid and the upper wall 
oscillate with the same frequency but with a finite phase difference , which is in 



agreement with the MD results by 



Thalakkottor and Mohsenil ( 2012 ) 



The continuum solution for oscillatory slip flows Eq. ([8]) was derived assuming 
constant slip lengths at the upper and lower walls. However, when analyzing the 
velocity profiles computed from MD simulations, we noticed that the fitting pa- 
rameters Li and L2 in Eq. ([8]) depend on the interfacial shear rate. The variation 
of the slip length as a function of the local shear rate computed at the upper 
and lower walls is plotted in Figs. [8] and [9] for different oscillation frequencies 10, 
amplitudes U, and times u>t n - For comparison, the data for steady shear flows are 
also presented in Figs.[8]and[9]on the log-linear scale to emphasize the low shear 
rate region. In all cases, the slip length for both oscillatory and steady-state flows 
is nearly constant at low shear rates jr < 0.01 and it increases linearly (see inset 
in Fig. [5} at higher shear rates. 

The deviation from the steady-state results in Fig. [8] is most pronounced at 
the highest frequency ujt = 0.1 and the largest amplitude U = 2.0 cr/r, when 
the magnitude of the wall acceleration is maximum, i.e., when ut = and n [see 
Tableland Fig.[5](a)]. Notice also that the MD velocity profiles near the upper 
wall develop pronounced oscillations at ut = 7r/2, and Zir/A in Fig.[5](a). It 
can be further observed that, the data in Fig.|8](a) are scattered at low shear rates 
because the velocity profiles near the lower wall are not significantly affected by 
the oscillating upper wall at the highest frequency ut = 0.1 (see Fig. [3]), and the 
statistical errors due to thermal fluctuations become relatively large. Similarly, the 
averaged fluid velocity is nearly zero when ut = and tt at the lowest frequency 
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LOT = 10" 4 (see Fig.ED, and, as a result, both L s and 7 are subject to statistical 
uncertainty when 7T < 0.005 in Fig. 0(b). Remember that, during each cycle, the 
data were averaged for the time interval T/100, and, thus, significantly longer 
averaging time is required to resolve accurately the velocity profiles in oscillating 
flows. It is expected, however, that with further averaging, the data in Figs. [8] and 
[9] for oscillatory flows at low shear rates will converge to the steady-state results. 

In the case of slip flow over a planar, impermeable solid surface, the friction 
coefficient that relates the wall shear stress and slip velocity is equal k = fi/L s , 
when the slip le ngth is computed by linear extrapolation of the velocity profile 



to zero velocity (jWillmott and Tallon 



2007). However, at higher frequencies, as 



shown in Fig. [51 the slip lengths L\ and L2 computed using Eq. ([8J deviate from 
the slip length in steady-state flows, and, therefore, these values do not provide 
an accurate estimate of the friction coefficient. Also, the estimate of the interfacial 
shear rate and the corresponding slip length directly from the MD velocity profiles 
(e.g., Figs.[3j[7|) is not precise because of the slight nonlinearity of the velocity 
profiles near interfaces and the ambiguity in choosing the size of the fitting region. 
To avoid the uncertainty associated with fitting the velocity profiles, the friction 
coefficient in oscillatory flows was estimated from the relation a xz (t n ) = ku s (t n ). 
The wall shear stress a X z(t n ) was computed as a ratio of the total tangential force 
between the fluid monomers and wall atoms to the wall area, and then averaged 
over the time interval T/100. At the same time, the slip velocity was calculated 
from the velocity and density profiles as follows: 



rz2 , rz2 

Us = u x (z)p(z)dz I I p(z)dz 

J z 1 J z 1 



(12) 
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where the integrals were taken over the width of the first peaks in the density pro- 
files. Naturally, the fluid slip velocity at the oscillating upper wall is the difference 
between the velocity of the adjacent fluid layer and the upper wall speed. 

As was shown in Fig.[2](a), in the presence of a solid substrate, the fluid 
monomers tend to form several distinct layers that are gradually decaying to a 
uniform bulk density. In addition to the density layering, the periodic surface po- 
tential typically induces an in-plane order within the adja cent fluid layers provided 



that the wall-fluid interaction energy is sufficiently high (jThompson and Robbins 



19901) . The characteristic signature of such ordering is the appearance of several 
sharp peaks in the static structure factor, which is defined as follows: 

E 

0=1 



1 I 1 |2 

s ^ = w f \E"' Tj \ ' (13) 



where the sum is taken over Nt fluid monomers in the first layer and rj = (xj , yj) is 
the position vector of the fluid monomer. These peaks are most pronounced at the 



first r eciprocal lattice vectors of the underlying substrate (jThompson and Robbins 



1990). Examples of the averaged str ucture factor an d its dependence on the slip 



velocity were previously reported by 



Priezievl (|2007al) for a similar computational 



setup. More recently, it was shown for several liquid-on-solid systems that the 
friction coefficient in steady flows correlates well with the product of the no rmalized 
peak in the structure factor and the contact density of the first fluid layer (jPrieziev 



2010|). 

In the present study, the friction coefficient in oscillatory and steady-state flows 
is plotted in Figs.llOland llll as a function of the combined variable S(0)/ [S(G\) p c ], 
where Gi = (9.04 er -1 , 0) is the first reciprocal lattice vector in the flow direction. 
For both types of flows, the friction coefficient and the induced fluid structure are 
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reduced at larger slip velocities. As is evident, the friction coefficient extracted 
from oscillatory flows agrees well with the the steady-state values, except that the 
data for periodic flows are scattered as small slip velocities, which is similar to 
the rate dependence of the slip length reported in Figs.[8]and[9] Interestingly, the 
agreement at higher frequencies is much better for the friction coefficient (shown 
in Fig.llOp than for the slip length in Fig.[8] Note also that at the highest frequency 
lot = 0.1, the period of oscillation T = 62.83t is about two orders of magnitude 
larger than the typical oscillation time of the LJ monomers; but nevertheless, 
the structure factor and the contact density are nearly the same as in steady- 
state flows. These results suggest that slip boundary conditions for high-frequency 
oscillatory flows are more accurately described by the dynamic friction coefficient 
rather than the slip length as a function of shear rate. 



4 Conclusions 

In this paper, we have investigated steady and oscillatory Couette flows with 
slip boundary conditions using molecular dynamics simulations. In both cases, 
the laminar flows were induced by the moving upper wall while the lower wall 
remained stationary. The simulations were performed in a wide range of oscillation 
frequencies; namely, when the Stokes boundary layer thickness is smaller than 
the channel width at the highest frequency, and, on the other hand, at lower 
frequencies that correspond to quasi-steady flows. For the chosen parameters, the 
liquid and solid phases form incommensurate structures at the interface, which is 
characterized by a finite slip length that increases almost linearly with shear rate. 



16 



Nikolai V. Priczjcv 



We found that the velocity profiles computed in MD simulations are well de- 
scribed by the Stokes flow solution with the slip length as a fitting parameter that 
depends on the local shear rate. The rate dependence of the slip length obtained 
in steady-state shear flows is reproduced in oscillatory flows when the slip length 
is measured as a function of the absolute value of the local shear rate. The MD 
data for oscillatory flows at low shear rates are relatively noisy, which, however, is 
not surprising given that the velocity profiles were averaged over a small fraction 
of the period during each cycle. For both types of flows, the friction coefficient at 
the liquid-solid interface correlates well with the structure factor and the contact 
density of the first fluid layer. 

Acknowledgements Financial support from the National Science Foundation (CBET-1033662) 
is gratefully acknowledged. Computational work in support of this research was performed at 
Michigan State University's High Performance Computing Facility. 



References 

Allen MP, Tildesley DJ (1987) Computer Simulation of Liquids. Oxford University Press, New 
York 

Asproulis N, Drikakis D (2010) Boundary slip dependency on surface stiffness. Phys Rev E 
81:061503 

Barrat JL, Bocquct L (1999) Influence of wetting properties on hydrodynamic boundary con- 
ditions at a fluid/solid interface. Faraday Discuss 112:119-127 

Bocquct L, Barrat JL (2007) Flow boundary conditions from nano- to micro-scales. Soft Matter 
3:685-693 

Du B, Goubaidoullinc I, Johannsmann D (2004) Effects of laterally heterogeneous slip on the 
resonance properties of quartz crystals immersed in liquids. Langmuir 20:10617-10624 



Oscillatory Couette flows with slip boundary conditions 



17 



Ellis JS, Hayward GL (2003) Intcrfacial slip on a transverse-shear mode acoustic wave device. 

J Appl Phys 94:7856-7867 
Ferrante F, Kipling AL, Thompson M (1994) Molecular slip at the solid-liquid interface of an 

acoustic-wave sensor. J Appl Phys 76:3448-3462 
Hansen JS, Ottesen JT (2006) Molecular dynamics simulations of oscillatory flows in microflu- 

idic channels. Microfluid Nanofluid 2:301-307 
Hansen JS, Daivis PJ, Todd BD (2007) Local linear viscoelasticity of confined fluids. J Chcm 

Phys 126:144706 

Kannam SK, Todd BD, Hansen JS, Daivis PJ (2012) Slip length of water on graphenc: Limi- 
tations of non-equilibrium molecular dynamics simulations. J Chem Phys 136:024705 

Karniadakis GE, Beskok A, Aluru N (2005) Microflows and Nanoflows: Fundamentals and 
Simulation. Springer, New York 

Khaled ARA, Vafai K (2004) The effect of the slip condition on Stokes and Couette flows due 
to an oscillating wall: exact solutions. Int J Nonlinear Mcch 39:795-809 

Kharc R, de Pablo J, Yethiraj A (2001) Molecular simulation and continuum mechanics in- 
vestigation of viscoelastic properties of fluids confined to molecular ly thin films. J Chcm 
Phys 114:7593-7601 

Li Y, Xu J, Li D (2010) Molecular dynamics simulation of nanoscale liquid flows. Microfluid 
Nanofluid 9:1011-1031 

Matthews MT, Hill JM (2009) On three simple experiments to determine slip lengths. Mi- 
crofluid Nanofluid 6:611-619 

Ng CO, Wang CY (2011) Oscillatory flow through a channel with stick-slip walls: complex 
Navicr's slip length. J Fluid Eng 133:014502 

Niavarani A, Priezjev NV (2008a) Rhcological study of polymer flow past rough surfaces with 
slip boundary conditions. J Chem Phys 12:144902 

Niavarani A, Priezjev NV (2008b) Slip boundary conditions for shear flow of polymer melts 
past atomically flat surfaces. Phys Rev E 77:041606 

Niavarani A, Priezjev NV (2010) Modeling the combined effect of surface roughness and shear 
rate on slip flow of simple fluids. Phys Rev E 81:011606 



18 



Nikolai V. Priezjev 



Pahlavan AA, Freund JB (2011) Effect of solid properties on slip at a fluid-solid interface. 
Phys Rev E 83:021602 

Priezjev NV (2007a) Rate-dependent slip boundary conditions for simple fluids. Phys Rev E 
75:051605 

Priezjev NV (2007b) Effect of surface roughness on rate-dependent slip in simple fluids. J 
Chem Phys 127:144708 

Priezjev NV (2009) Shear rate threshold for the boundary slip in dense polymer films. Phys 
Rev E 80:031608 

Priezjev NV (2010) Relationship between induced fluid structure and boundary slip in 
nanoscale polymer films. Phys Rev E 82:051603 

Priezjev NV (2011) Molecular diffusion and slip boundary conditions at smooth surfaces with 
periodic and random nanoscale textures. J Chem Phys 135:204704 

Priezjev NV (2012) Interfacial friction between scmiflcxiblc polymers and crystalline surfaces. 
J Chem Phys 136:224702 

Priezjev NV, Troian SM (2004) Molecular origin and dynamic behavior of slip in sheared 
polymer films. Phys Rev Lett 92:018302 

Priezjev NV, Darhuber AA, Troian SM (2005) Slip behavior in liquid films on surfaces of pat- 
terned wettability: Comparison between continuum and molecular dynamics simulations. 
Phys Rev E 71:041608 

Priezjev NV, Troian SM (2006) Influence of periodic wall roughness on the slip behaviour at 
liquid/solid interfaces: molecular-scale simulations versus continuum predictions. J Fluid 
Mech 554:25-46 

Qian TZ, Wang XP, Sheng P (2005) Hydrodynamic slip boundary condition at chemically pat- 
terned surfaces: A continuum deduction from molecular dynamics. Phys Rev E 72:022501 

Thalakkottor JJ, Mohseni K (2012) Analysis of boundary slip in a flow with an oscillating 
wall. larXiv:1207.7090l 

Thompson PA, Robbins MO (1990) Shear flow near solids: Epitaxial order and flow boundary 

conditions. Phys Rev A 41:6830-6837 
Thompson PA, Troian SM (1997) A general boundary condition for liquid flow at solid surfaces. 

Nature 389:360-362 



Oscillatory Couette flows with slip boundary conditions 



19 



Wang FC, Zhao YP (2011a) Slip boundary conditions based on molecular kinetic theory: The 
critical shear stress and the energy dissipation at the liquid-solid interface. Soft Matter 
7:8628-8634 

Wang FC, Zhao YP (2011b) The unique properties of the solid-like confined liquid films: A large 
scale molecular dynamics simulation approach. Acta Mcchanica Solida Sinica 24:101-116 

Willmott GR, Tallon JL (2007) Measurement of Newtonian fluid slip using a torsional ultra- 
sonic oscillator. Phys Rev E 76:066306 

Yang SC, Fang LB (2005) Effect of surface roughness on slip flows in hydrophobic and hy- 
drophilic microchannels by molecular dynamics simulation. Molecular Simulation 31:971- 
977 

Zhang HW, Zhang ZQ, Ye HF (2012) Molecular dynamics-based prediction of boundary slip 
of fluids in nanochannels. Microfluid Nanofluid 12:107-115 



20 



Nikolai V. Priczjev 



Table 1 The oscillation frequency a) (in units r~ ), the oscillation period T = 2tt/u) (in units 
r), the maximum amplitude of the upper wall velocity (in units ct/t), the upper estimate of the 
Reynolds number Re max = AUhp/fi, the Stokes boundary layer thickness 5 = fi/poJ (in 
units <t), and the corresponding Reynolds number Re max = Af7<5p//x, when 8 < h = 23.58 it. 
In the definition of the Reynolds numbers, AC/ is the maximum variation of the tangential 
fluid velocity component across the channel. 





T/r 


U max 




S/(T 


Re s 

1 max 


0.1 


62.83 


2.0 


6.5 


7.29 


2.0 


0.01 


628.32 


4.0 


18.8 


23.04 


18.4 


0.001 


6283.19 


6.0 


27.6 


72.86 




0.0001 


62831.85 


(5.0 


27.6 


230.4 





Table 2 The instantaneous slip lengths and shear rate magnitudes at the oscillating upper 
wall obtained from the best fit of the MD velocity profiles in Fig.[5]using Eq. 0. The amplitude 
of the upper wall velocity is U = 2.0 ct/t for u)t = 0.1 and U = 4.0<t/t for ujt = 0.01. The 
same data as in Fig. [8] 

ujt 7r/4 tt/2 3tt/4 t 5tt/4 3tt/2 7tt/4 

ujt = 0.1 

7t 0.037 0.124 0.136 0.071 0.037 0.125 0.136 0.071 

L s /a 9.8 9.8 10.1 10.0 9.8 9.7 10.1 10.0 

ujt = 0.01 

7t 0.080 0.140 0.129 0.045 0.080 0.140 0.129 0.045 

L s /a 7.6 10.3 11.1 9.5 7.6 10.3 11.1 9.5 
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Fig. 1 (Color online) Positions of fluid monomers (open blue circles) and wall atoms (filled 
gray circles). The upper wall oscillates with the angular frequency u> in the x direction (indi- 
cated by the double-sided arrow), while the lower wall is always stationary. 




z/o 

Fig. 2 (Color online) Ensemble-averaged (a) density and (b) velocity profiles for the indicated 
upper wall speeds and steady-state flow conditions. The vertical dashed lines at z = ±11.79 a 
indicate the reference planes for computing the slip length. The vertical axes at z = ±12.29 <r 
coincide with the location of the fee lattice planes which are in contact with the fluid phase. 
The inset shows the slip length as a function of shear rate when u) = 0. 
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Fig. 3 (Color online) Averaged velocity profiles for oscillatory flows with amplitudes (a) U = 
0.25 cr/r and (b) U = I.Oct/t and frequency lot = 0.1 at times uit = 0, 7r/4, n/2, 37r/4, and7r. 
The red dashed curves are the least square fits of Eq. {5J to the MD data. The vertical dashed 
lines at z = ±11.79 cr indicate the reference planes for computing the slip length and shear 
rate from the best fits to Eq. JS). 
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Fig. 4 (Color online) Velocity profiles for time-periodic flows with amplitudes (a) U = 
0.25 cr/r and (b) U = 3.0<r/r and frequency lot = 0.01 at times uit = 0, tt/4, tt/2, 3 tt/4, and7r. 
The red dashed curves are the best fits of the MD data using Eq. {5J . The vertical dashed lines 
at z = ±11.79(7 denote the location of liquid-solid interfaces. 
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Fig. 5 (Color online) Averaged velocity profiles for (a) u)t = 0.1 and U = 2.0cr/r and (b) 
lot = 0.01 and U = A.Oct/t at times tot = 0, 7r/4, 7r/2, 3 7r/4, and-n - . The dashed curves are 
the best fits of the MD data using Eq. JS}. The corresponding slip lengths and shear rate 
magnitudes at the oscillating upper wall are listed in Table|2] 
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Fig. 6 (Color online) Velocity profiles for oscillatory flows with amplitudes (a) U = 0.25 it/t 
and (b) U = 5.0<t/t and frequency lot = 0.001 at times ut = 0, 7r/4, tt/2, 3tt/4, and7r. The 
dashed curves represent the best fits of the MD data using Eq. (5). The vertical axes at z = 
±12.29cr coincide with the location of the fee lattice planes. 




Fig. 7 (Color online) Velocity profiles for time-periodic flows with amplitudes (a) 
U = 0.25 a It and (b) U = 6.0ct/t and frequency lot = 0.0001 at times tot = 
0, 7r/4, 7r/2, 3 7r/4, and7r. The red dashed curves indicate the best fits of the MD data using 
Eq. ©. 
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Fig. 8 (Color online) The slip length L a /a as a function of shear rate for (a) lot = 0.1 and 
(b) lot = 0.01. The slip length and the interfacial shear rate are evaluated at the stationary 
lower wall (o) and at the oscillating upper wall (o). The data for the largest amplitudes U 
are presented in Table[2] The data for steady shear flows (□) are the same as in the inset of 
Fig. [2](a). The black curves are guides to the eye. 
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Fig. 9 (Color online) Variation of the slip length as a function of shear rate for (a) u)t = 0.001 
and (b) lot = 0.0001. The slip length and shear rate are computed at the stationary lower wall 
(o) and at the oscillating upper wall (o). The data (□) are the same as in the inset of Fig.[2](a). 



,_, 6 



CO 



CO 



JO 



,_, 6 



(a) 



(b) 



32 




cox = 0.1 




cox = 0.01 



46 



60 



S(0) [StG^pol" 



Fig. 10 (Color online) Log-log plot of the inverse friction coefficient as a function of 
5(0)/ [iS(Gi) p c ] for (a) ojt = 0.1 and (b) lot = 0.01. The data for oscillatory flows are denoted 
by the red diamonds (oscillating upper wall) and blue circles (stationary lower wall) . The data 
for steady shear flows are indicated by the black squares. The black curves are guides to the 
eye. 
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Fig. 11 (Color online) The inverse friction coefficient versus S(0)/ [S(G\) p c ] for (a) lot = 
0.001 and (b) lot = 0.0001. The friction coefficient and fluid structure are computed at the 
oscillating upper wall (o) and at the lower stationary wall (o). The data for steady shear flows 
(□) are the same as in Fig. 1101 



